stw <- function (ref, samp = NULL, selected.traces,
		 time = 1:ncol(ref)-0.5,
		 bspl = 40, spldeg = 3, 
		 init.coef = rep(0, bspl),
		 lambda = 1, kappa = 0.01, 
		 diffdeg = 2,  try = FALSE,
		 warp.type = c("individual", "global"),
		 smooth.param = ifelse(try, 0, 1e05),
                 verbose = FALSE, ... )
{
  if(length(init.coef) != bspl) {
    warning(paste("Number of initial coefficients does not match number of B-splines; setting 'bspl' to ", 
		  length(init.coef), sep=""), immediate = TRUE)

    bspl <- length(init.coef)
  }

  stw.par.list <- list(bspl = bspl, spldeg = spldeg, lambda = lambda, kappa = kappa, diffdeg = diffdeg)

  if (any(class(ref) == "ptw")) {
    ptwList <- ref
    ref <- ptwList$reference
    samp <- ptwList$sample
    time <- ptwList$warp.fun
    ptw.w.deriv <- ptwList$deriv.warp.fun
    warp.type <- ptwList$warp.type
  } else {
    ptw.w.deriv = 1 # No ptw derivative obviously, but simply the derivative of w(t) = t
    warp.type <- match.arg(warp.type)
  }

  if (is.vector(ref)) ref <- matrix(ref, nrow = 1)
  if (is.vector(samp)) samp <- matrix(samp, nrow = 1)
  if (is.vector(time)) time <- matrix(time, nrow = 1)
  if (nrow(ref) > 1 && nrow(ref) != nrow(samp))
    stop("The number of references does not equal the number of samples")
  if (nrow(samp) == 1) warp.type <- "individual"
  
  r <- nrow(samp)
  
  if (!missing(selected.traces)) {
    samp <- samp[selected.traces,, drop=FALSE]  

    if (nrow(ref) > 1)
      ref <- ref[selected.traces,, drop=FALSE]
          
  }

  if (is.vector(init.coef)) init.coef <- matrix(init.coef, nrow = 1)
  if (warp.type == "global") {
    if (nrow(init.coef) != 1)
      stop("Only one warping function is allowed with global alignment.")
  } else {
    if (nrow(init.coef) != nrow(samp))
      if (nrow(init.coef) == 1) {
        init.coef <- matrix(init.coef, byrow = TRUE,
                         nrow = nrow(samp), ncol = length(init.coef))
      } else {
        stop("The number of warping functions does not match the number of samples")
      }
  }

  if (is.vector(time)) time <- matrix(time, nrow = 1)
  if (warp.type == "global") {
    if (nrow(time) != 1)
      stop("Only one warping function is allowed with global alignment.")
  } else {
    if (nrow(init.coef) != nrow(samp))
      if (nrow(init.coef) == 1) {
        init.coef <- matrix(init.coef, byrow = TRUE,
                         nrow = nrow(samp), ncol = length(init.coef))
      } else {
        stop("The number of warping functions does not match the number of samples")
      }
  }

  if (warp.type == "individual") {
    w.deriv <- w <- matrix(0, nrow(samp), ncol(ref))   
    a <- matrix(0, nrow(samp), ncol(init.coef))
    v <- rep(0, nrow(samp))
    warped.sample <- matrix(NA, nrow=nrow(samp), ncol=ncol(samp))

    for (i in 1:nrow(samp)) {
      if (verbose & nrow(samp) > 1)
        cat(ifelse(nrow(ref) == 1,
                   paste("Warping sample", i, "with the reference \n"),
                   paste("Warping sample", i, "with reference \n", i)))

      if (nrow(ref) == 1) {
        rfrnc <- ref
      } else {
        rfrnc <- ref[i, , drop = FALSE]
      }

      if (nrow(time) == 1) {
        t <- time
      } else {
        t <- time[i, , drop = FALSE]
      }


      warp.res <- spmwarp(rfrnc, samp[i, , drop = FALSE], 
			  init.coef[i,], stw.par.list, 
			  t, try = try, 
			  smooth.param = smooth.param, ...)
      
      w[i, ] <- warp.res$w
      a[i, ] <- warp.res$a
      v[i] <- warp.res$v
      if (identical(dim(ptw.w.deriv), dim(w.deriv))) {
        w.deriv[i, ] <- warp.res$stw.w.deriv + ptw.w.deriv[i, ]
      } else {
	  w.deriv[i, ] <- warp.res$stw.w.deriv + ptw.w.deriv
      }
      warped.sample[i, ] <- interpol(w[i, ], samp[i, ]) * w.deriv[i, ]
    }
  } else {
    warped.sample <- matrix(NA, nrow=nrow(samp), ncol=ncol(samp))
    
    if (nrow(ref)==1) 
      ref <- matrix(ref, nrow = nrow(samp), 
      			ncol = ncol(ref), byrow = TRUE)
    
    if (verbose) {
      if (nrow(ref) == 1) {
        cat("Simultaneous warping of samples with reference... \n")
      } else {
        cat("Simultaneous warping of samples with references... \n")
      }
    }

    warp.res <- spmwarp.global(ref, samp, c(init.coef), stw.par.list, 
				 time, try = try, smooth.param = smooth.param,  ...)

    w <- as.matrix(warp.res$w)
    a <- t(as.matrix(warp.res$a))
    v <- warp.res$v
    stw.w.deriv <- as.matrix(warp.res$stw.w.deriv)
    w.deriv <- stw.w.deriv + ptw.w.deriv

    if (dim(w.deriv)[1] == 1) {
      warped.sample <- t(sapply(1:nrow(samp),
			  function(i) {
			    interpol(w, samp[i, ]) * w.deriv
			  }))
    } else {
	warped.sample <- t(sapply(1:nrow(samp),
			  function(i) {
			    interpol(w, samp[i, ]) * w.deriv[i, ]
			  }))
    }
  }
  
  if (verbose) cat("\nFinished.\n")  

  result <- list(reference = ref, sample = samp,
                 warped.sample = warped.sample,
                 warp.coef = a, warp.fun = w,
                 crit.value = v, warp.type = warp.type, 
		 deriv.warp.fun = w.deriv)
  class(result) <- c("stw", "ptw")
  result
}